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ABSTRACT; 


General  formulas  for  the  first-order  effects  of  the 
oblateness  of  the  Earth  upon  the  orbits  of  satellites 
were  reported  in  1957  by  R.   E.   Roberson  and  D.    G. 
King-Hele.      Their  derivations  employed  methods  which 
are  unfamiliar  to  those  not  having  studied  the  mathe- 
matics of  astronomy.     In  this  Technical  Report  / 
Research  Paper  No.    76  s   an  elementary  derivation  is 
given  for  these  useful  and  important  results. 


PRESENTATION 


This  paper  will  be  presented  26  May  1967  at  the 
Canadian  Congress  of  Applied  Mechanics  s   Quebec. 


An  Elementary  Derivation  of  the  Principal  Effects   of  the 
Oblateness   of   the   Earth   upon   the   Orbits   of  Satellites 

by 

John  E.    Brock 

Professor  of  Mechanical  Engineering 

Naval   Postgraduate   School 

Monterey,      California 

Independently  of  each   other,    in   1957  R.    E.    Roberson  in   the 
United   States    and  D.    King-Hele   in  Great  Britain  provided   general 
formulas   giving   the    first   order  effects   of  the   oblateness   of  the 
Earth   upon  satellite   orbits,   previous   analyses   by   others  having 
contained   inconvenient   restrictions.      Roberson   and  King-Hele   obtained 
their  results  by   the   use   of  methods  which   are   familiar  to   those  who 
are   trained   in   the  mathematics   of   Astronomy  but  not    familiar   to 
those  having  onlv   the   usual   training   in  Engineering  Mechanics,  who, 
accordingly,    do  not  have  easy   access   to   these   useful   and   important 
results . 

In   the  present  paper   these    results    are   derived  by   a  quite  ele- 
mentary procedure   requiring  nothing  more    than  elementary  geometry, 
vector  algebra,    and  mechanics.      An   appendix  considers   a   few  simpli- 
fied models   of   the   Earth  which   have   the   same    first   order  effects   as 
does    an      oblate   Earth.      This  material   should  be   of   interest    to  en- 
gineers  engaged   in    the   space   effort   as    a  means   of  understanding 
the  primary  perturbations  which   result   from  oblateness. 

A  complete    table   6f  notation   is   given   in  Appendix  B.      This 
table   is  not   arranged   alphabetically,  however,   since   so  many  of  the 
definitions   of  the   quantities   depend   on   others    in    the    list.    Instead 
the   arrangement  is   in   "operational   sequence",    that  is,    in   the   order 


in  which   the   symbols    logically  occur  in   the   development.      The 
notation   agrees   generallywith   that  used   in   astronautical   litera- 
ture;   for  example,   V  means   geopotentlal ,   v  means    true   anomaly,    and 
V  denotes   satellite    (vector)    velocity. 

First  we   establish   some   coordinate   systems.    0  is    the  mass 
center  of   the  Earth   and  n   is    a  unit   vector  parallel    to   the   axis   of 
rotation  of  the  Earth  pointing  North,      a   and  b   are   fixed   unit   vec- 
tors  in  the  Earth's   equatorial  plane   so  as    to   form  a  right  handed 
cartesian   system  a,   b,   n;    actually,    of   course,  because   of  the  pre- 
cession  of   the   equinoxes    and   other  slow  motions,    these   vectors    are 
not   truly   fixed,   but  we   may   regard   them  as  being   fixed   for   the  pur- 
poses  of   this    analysis.      We   could   take   a   to  be   in   the   direction  of 
the  projection  upon   the  Earth's   equatorial  plane   of  the   line   from 
0  to  the   first  point  of  Aires,  but   there   seems    to  be  no  uniformity 
in    this  matter.      We   observe    that   conventionally   the   orbits   of  Earth 
satellites    are    referred   to   the  Earth's   equatorial  plane   in   the   same 
way   that  planetary   orbits   are   referred   to   the  plane   of  the  ecliptic, 
but   there   does   not   seem  to  be   perfect   uniformity  with   regard   to   a 
reference  position    for  longitude   such   as    that  provided  by   the    first 
point   of  Aires. 

Let  Q  be   the  instantaneous  position  of  the   satellite,  with   r  ■ 
OQ  and   V  =   r,   the   dot  denoting  differentiation  with   respect   to  time. 

We   define   the   scalar  coordinates   of  satellite  position  by   the  equations 

x  =  r«a,       y  ■  r»b,       z  ■  r«n  [la,b,c] 

so  that 

r  -  xa  +  yb  +  zr\  [2] 


The  geopotential   for  an   oblate   Earth    (see   Appendix  A)   may  be 
approximated  by   the  expression 

V  -  -ur~l[l  +  Jr"2(l  -   3s2r~2)/3]  [3] 

where  y   ■  CM     is    the  product  of   the   universal    constant   of  gravitation 
by   the  mass   of   the   Earth,   J  is    the   coefficient   of  the   second   zonal 


harmonic,    and  r  =   /r»r.      An   approximate  numerical   value    for  \i    is 

9.66*10u    (statute  miles)    (seconds)      .      At   the   surface   of   the  Earth, 

_2 
at   the  equator  where   r  =   3964   statute   miles,    the   product  Jr    '  has 

-3 
the   value    1. 623*10        (diraensionless) .      Thus,   evidently,    for  an  Earth 

satellite,  the  second  term  in  the  square  brackets  in  Eq .  [3]  is  much 
smaller  than  unity.  For  a  further  discussion  of  Eq.  [3],  see  Appen- 
dix A. 

The    force,   per  unit  mass   of  satellite,   exerted  on   the   satellite 
by   the  gravitational   attraction  of  the   Earth   is 

*  i  t,  JUL     k^V        3V  r/1 

f  =  -  grad  V  =  -  a^  -b-  -n-  [4] 

2  2  2  2 

and  to  evaluate   this   expression,  we  note    first   that  r     =  x     +  y     +  z 

so    tnat  n  o        /o  o 

.   n      ,  n        2       ,,   2Nn/2    .    2 
lL  _  QL  -111  _  d(r  )        .8r 
3x  *  dr»2   dx  "  dr2        3x 

=    (n/2)(7»2)n/2_1(2x)   =  nxrn~2  [5] 

with  similar  formulas   in  which  y   and  a   replace  x.     Thus 
f  =    ugradfr"1  +  JJ*~3/3  -  J32r~5] 

=  v  [(-xr~3  -  Jxr"5  +  5J^32r~7)a  + 

(-yr      -  Jyr      +  5Jyz  r     )b  + 

-3  -5  2  -7  -5 

(-zr       -  Jzr>      +  5Jzz  r       -  23zv     )n] 

and  by  using  Eq.    [2]    this  may  be  more  briefly  written   as 


f  -  -  vr~3r  -  Jur~5r(l  -  5z2r>~2)  -  2Jyzr"5n 

_2 
-  -  yr     e     -  AJue     -  BJun  [6] 

where    *    -   r/r  is   a  unit  vector  parallel   to   r  and 
r 

A  -  r~A(l  -  532r~2),        B  =  2sr~5  [7a, b] 

Note  that  if  J  were  zero  in  Eq.  [7] ,  we  would  have  the  simple  inverse 

square  law  of  attraction  toward  the  mass  center  of  the  Earth. 

We  now  define  h  to  be  the  angular  momentum  (with  respect  to  0) 

of  a  unit  mass  of  satellite,  i.  e., 

h  -  r*v  [8] 

Since 

v  =  f  -  d(re  )/dt  ■  re  +  re 
r        r     r 

and  since  e   *e  ■  0,  we  have 
r     r 

2 
r*h  =  rx(rxv)  -  r«vr  -  r>rv  =  r  (e   »ve     -v) 

r       r 

■  r  [(!•  +  0)e     -  re     -  re   ]  =  -r  e  [9] 

r  r  r  r 

We  now  use  the  fundamental  principle  of  Mechanics  that  the  rate  of 
change  of  angular  momentum  (about  0  which  we  consider  to  be  a  fixed 
point  since  we  are  concerned  with  the  motion  of  0  relative  to  0)  is 
equal  to  the  moment  about  0  of  the  applied  force;  that  is 

fi  -  rxf  [10] 

In  the  present  case  with  only  one  moving  particle,  this  formula  is 
very  easy  to  obtain  from  Eq.  [8] ,  from  which 

h  «  f*v  +  r*v  -  r*v  =  rxf 
since  V  is  the  acceleration  which  is  equal  to  f ,  the  force  per  unit 
mass. 

Next,  consider  the  expression  V*h.   Taking  the  derivative  with 
respect  to  time  and  using  the  notation  e,  ■  h/7z,  where  h   «  A)*\) ,  we 

have 

d 

~(vxh)  -  vxh  +  vxh  -  fxh  +  vx(rxf) 


_3 
■  -  pr     r*h  -  Jp(Ae    +  Bn)xh  +  vx[rx(-BJun)] 

-l 
-  ye    +  uJMAe,  *e    +  Be  xn  +  Bfc     vx(nxr)]  [111 

r  h      r  h 

Integrating   this   expression  with   respect    to   time   gives 

Vxh   =  u(e     +  S)  [12] 

where   S   is    a  vector  whose    time-derivative   satisfies 

•  -1 

s  ■  jft[Ae.  *e    +  Be,  *n  +  Bfr     vx(nxr)]  [13] 

h     r  n 

Also,    from  Eq.     [10]    and   Eq.     [6], 

h  =  BJynxr  [14] 

Since  we  have   seen,    in  Eq.    [3],    that   the   terms    involving  J  are   small 
compared  to   those  not   containing  J,  we   regard   J  as   a   "small"   quantity. 
Thus,   Eq.    [13]    and  Eq .    [14]    indicate   that   S    and  h   are    "nearly"   con- 
stant vectors;    that   is,   their  magnitudes   and   orientations   are  but 
slowly   changing. 

Next,    consider  the   scalar  quantity 

2 
r»v*h  ■  rxvh  =  h    ■  u(r  +r*s) 

=  ur(l  +  e   »S)  -  ur(l  +  s   costf)  [15] 

where   s  =   /§»S    and   cos   V  =  e    'S/s.      V   is   the   angle,   shown   in  Fig.    1, 

which   is  known   as    the    true   anomaly.      If  we   define   the  quantity 

p  =  h2M~l  [16] 

which  has  the  dimensions  of  a  length,  we  may  write 

p   »  r(l  +  s  cos  v)  [17] 

-1*     • 
In  Eq.  [17]  only  r>   and  V   vary  "rapidly",  p  =  2h]i     h   and  s  both 

involving  J  as  a  factor.   Eq.  [15]  is  that  of  an  ellipse,  called  the 

instantaneous  elliptical  orbit  (IEO) ,  from  which  the  actual  orbit  of 

the  satellite  begins  to  differ  as  time  goes  on.  The  point  on  this  IEO 

(point  P  in  Fig.  1)  corresponding  to  v   ■  0  is  called  the  instantaneous 


perigee  point,  toward  which  the  unit  vector  e  ■  S/e  points.   Also, 

the  eccentricity,  s,  of  the  IEO  is  called  the  instantaneous  eccentricity, 

The  angle  i   (0  ■  i   *w)  ,  satisfying  cos  i   ■  n»e,  ,  is  the  (instantaneous) 

inclination  of  the  orbit  plane.   A  unit  vector  in    the  direction  e  *n , 

i.  e. ,  the  vector  j  «  ev*n  csc  ^»  being  perpendicular  to  both  the  plane 

of  the  IEO  and  the  equatorial  plane,  thus  points  toward  the  descending 

node  D  (Fig.  1)  of  the  former;  -j  points  toward  the  ascending  node  A. 

The  angle  u>  *  <A0P  is  called  the  argument  of  perigee.   The  quantity  p 

is  seen  to  be  the  semi-latus  rectum  of  the  IEO. 

We  will  use  the  unit  vector  i  ■  j*e,  such  that  i,  j,  e,  form  a 

h  h 

right  handed   triad,   slowly   changing   in   orientation  since   j   and  e      are 

not   constant.      We   let  u  and  w  be   the   angles     BOQ  and     BOP,   respectively 

(Fig.    1)    so   that   cos   u  ■   l"  *e    ,    cos  w  ■   i  *e    .      Note   that  u  m  V  +  W  and 

r'  s 

that  w  is   "nearly"  constant  since   the   IEO  changes   only  slowly.      This 
implies   that  u  =  v,   approximately,   a  result  which  will  be   used  repeat- 
edly.     Note   also  that   u  =  w  +  tt/2  .      We   also  introduce   the   unit  vector 

eo  ■  e,  xe     so  that  e    ,   ert,  e,    also   form  a  slowly  moving  unit   triad. 
8hs  s        6       h 

Finally,  we   introduce   the  unit  vector  I   =  jxn   so  that   I ,  j ,  n    is   a  unit 

triad  which   rotates   slowly  about  n«      Thus,   in  particular,  we  have 

i    -  I    cos  i  +  n  sin  i  [18] 

e  *   i    cos  u  +  j   sin  u  [191 

r 

A  littie  manipulation  easily  establishes   the   following  relations, 
nxr  -     (-1  sin  u  +  j  cos  i  cos  u)  [20] 

z  -  re   »n  -  r  sin  i  cos  u  [21] 

r 

Now  it  is  obvious  that,  to  the  approximation  with  which  we  are  working, 

we  have 

e  -  we  *e  t22l 

r    h  r 

but  it  is  also  easily  possible  to  obtain  this  result  formally  as  follows 


e     ■   1    cos  u  +  j   sin  u  -    (1    sin  w  -  J   cos  u)u 

•  •  • 

-   (J   cos  u  -  1    sin  u)u 

since   1    and  j   nearly  vanish.      However, 

e.  «?     ■  e,  xj   sin  u  +  e,  x1    cos  u  ■  -i    sin  w  +  j   cos  u 
n     r         n  n 

and  Eq.    [22]    is   established.      From  this  we  have 

h  =.  rxf  =  rx(re    +  rue  xe  )  =  rue  x(e  *e  )  -  r  ue,_ 
r  n     r  r    '   n     r  h 

ft  =  r2w  [22a] 

r>~  =  r     h~  u  ■  p~  h~  u(l  +  e  cos  v)  [23] 

Substituting  Eqs.    [7b],    [20],    [21],    and    [23]    into  Eq.    [14],   expanding   cos   v 

■  cos(u  -  w) ,    and  multiplying  by   dt,  we   obtain 

dh  «  2Jyp     7z       sin  i    (1  +3  cos  u  cos  w     +  s  sin  u  sin  w) 

2 
x(-I   sin  u  cos  u  +  J    cos  i   cos  u)du  ,_., 

[24] 

It   is   elementary  but   tedious   to  obtain  the   indefinite   integral  of  this 
expression   and  since   the   result  of  such  an  exercise   does  not  seem  to  be 
of  particular  value,  we  will  perform  the   integration   for  one   complete 
passage   only,   as  u  increases  by   2tt,.    to  obtain  what  we   designate   as   Ah, 
the  net   change   in  h  per  passage.     To  do   this  most   easily,   note   that   for 
non-negative   integers  M,   N 

.   sin  x  cos  x  6x  =  0  [25] 

unless  both  M  and  N  are  even.      We  easily  obtain 

Ah   -   jj\mh~lp~l   sin   2i  [26] 

Since  j.e,    -  0,    the  variation,  Aft,   in  the  magnitude   of  h   is   of  higher  order 

in  J.     Thus  we   also  have 

—2  —1  —2 

Aeh  -  3Jvvh     p       sin  2£  -  jJwp       sin  2i  [27] 

using  Eq.    [16],   and   from  the  definition  of  j  we  have 

Aj  -  A(e,xn  cac  i)  -  Aehxn  esc  i  +  e^n  Acsc  i 


-2  -2 

Aj   -  2Jnp     j*n   cos   i  +  JAcsc  i  -  2Jirp     I    cos   i  [28] 

since  evidently  we  must  omit  the  term  In  J  because  j.j  ■  1.     We  may  also 

see  In  another  way  that   there  Is  no  j-component  since   the  variation  In t 

is   of  higher  order  In  J,   Ah  being  perpendicular  to  the  plane   in  which 

angle  i  is  measured.      We   observe   that   the   line   of  nodes   rotates    (in  a 

sense  opposite  to  that  of  the  satellite  itself  if  i   <  tt/2)  with  mean 

angular  velocity 

B*  -  -2HJTrp"2P"1cos   i  [29] 

where  P  is   the  orbital  period,  which  differs   only  by   a  term  of  order  J 
from  the  value 

P  »  2Tra3/2u~1/2  [30] 

(a  Is   the  semimajor  axis   of  ellipse)  which  would  obtain  if  J  were   zero. 
We  use   the  asterisk   in  this    formula  and  later  to  indicate   a  mean  value 
in  the  sense   that  we  have   indicated.      Thus  we   can   finally  obtain  the  very 
important   formula 

BN  «  -nc  cos  i  [31] 

where 

C  -   (g/r>E)l/2(J/rl)(rK/a)7/2(l  -  s2)"2  [321 

-2 
in  which  r»  represents   the  equatorial   radius   of  the  earth  and  g  ■  yr 

represents   the   acceleration  of  gravity   at  the  surface.      We  have  used   the 

following   relations 

2a  "  r(y-0)  +  r(u-*)   "  P/(1+s>  +  P/U-«>   "  2p/(l-s2) 

,.-2.-1        .  -7/2,,      2N-2   1/2        .  -7/2,,      2N-2   1/2     -1 
2Jto     P       »  3a         (l-«  )     y         »  Ja         (1-8  )     g       r         -  C 

E 

Returning  now  to  Eq.    [13],  we  wish  to  express  all  vectors  in  terms   of 

the    (slowly  moving)    triad  e    ,   e       «..     We  have 

son 


e  ■  e  cos  v   +  eA  sin  v:      n  ■  e,  cos  i  +   i  sin  i 
r    s  6  h 

1  =  e  cos  w  -  en   sin  w;     v  =  f  =  re  +  rue,  *e 
s         6  r      h  r 

It  is  now  a  matter  of  straightforward,  if  tedious,  substitution  into 
Eq.  [13].   Some  intermediate  evaluations  are  shown  below. 

e  xe  -  -e  sin  v  +  e„  cos  v 

h   r     s  6 

e,  *n  =  e,  *1  sin  i   -  sin  i    (e  sin  w  +  en   cos  w) 
h       h  s  8 


n*e  = 
r 


s  e         h 

sin  i   cos  w  -   sin  i   sin  u     cos  i 


cos  y 


sin  y 


n 


=  -e  cos  i  sin  y  +  e„  cos  t  cos  y  +  e,  sin  i   sin  w 
s  0  h 

2  • 
vx(nxr)  -  rre  x(nxe  )  +  r  w(e,  xe  )x(nxe  ) 
r     r        h  r      r 


=  rr(e  sin  i   sin  u  sin  y  -  e.  sin  i   sin  k  cos  v   +  e,  cos  t) 
s  0  h 

2» 
+  r  w  sin  t  sin  w  (e  cos  v   +  e.  sin  y) 

s         9 

2« 
From  Eq.  [22a]  we  have  r  u  "  h   and  from  Eq.  [17]  and  Eq.  [23]  we  get 

•    3  -1  •  -1 

r>r   =  r  p     sv  sin  v  =  hs  sin  V   (1  +  s   cos  v)      .   Thus  we  can  write 

vx(nxr)  =  (1  +  s   cos  y)  Msin  t  sin  u  [e  (s  +  cos  y)  +  e  sin  y] 

+  e,  s  sin  y  cos  i) 
n 

-A     -2*  2 

We  also  note  that  we  may  write  Jhr       =  Jp     v(l  +  s   cos  y)  by  using 


•       • 


Eq.  [23]  and  Eq.  [17]  and  the  fact  that  u   =  V   approximately.   Thus,  upon 
substituting  into  Eq.  [13],  we  find 


•  -2  2  2 

S  -  Jp      (1+scos   y){(l-5sin  t   cos  u)(l  +  ecos  v)  (-e  sin  v  + 

2  . 
eacos   y)   +  2(l+scos   y)    sin  \   cos  u   (e  slnw+  e   cos  w)  + 
6  So 

2 
2   sin  i  sin  w  cos   w   [e   (s +  cos   v)  +  e  sin  v]   + 

s  9 

ee,  sin  2t  sin  y  cos  u}  v 

h 

Multiplying  by  dt  and  expanding  functions  of  u   in  terms  of  V     and 
w,  w   being  regarded  as  constant,  gives  an  expression  for  dS  which 


10 


may  be   integrated   for  general  V.     However,    It   is    the  net   change  per  period 

which  is   of  greatest  interest,   and  upon  use   of  Eq.    [25],  we    find   that 

-2  2 

AS   =*  JsTrp      [(2-3  sin  i)e     -  e,     sin  2i   sin  w]  [33] 

Since  S*AS  ■  0,   the  variation  of  8   is   of  higher  order  in  J,    and  we  have 

-2  2 

Ae     =  Ji\p      [(2-3  sin  ^)e     -  e     sin   2^  sin  w]  [34] 

For  the   argument   of  perigee,  we  have 

-   cos   co  ■  j  #e     *  sin  W 
s 

-   cos(w+Aw)    ■   -    cos    oj  +  Aw  sin   to 

-  (J+Aj)-(e  +Ae  ) 

s        s 

-  j«e    +  j«Ae    +  Aj«e 


s 
s " 


inW  +  jTrp   l(2-3sin  ^)j-«g  +  2l»er  cos   i] 


-2  2 . 

=   sin  W  +  Jirp    '    cos  W   (4  -   5   sin  t) 


Since   sin  to  ■  cos  w,  we   clearly  have 


Aw  =  Jirp   2(4  -   5   sin2i)  [35] 

Dividing  by  P,    the    time   of  passage,  we  have 

w*  =    (C/2)(4  -   5   sin2   t)  [36] 

for  the  mean  rate   of  change   of  argument  of  perigee.      From  Eq.    [31]   we  see   that 

.* 

fl     »  -  C  cos   t 

is   the  mean  rate   of  change   of  the   longitude   of   the   ascending  node.    Combining 
this   result  with   that  of  Eq.    [36],  we   find  that   the  mean  rate  of  change   of 
the   longitude   of  perigee   is 

w*  -    (C/2)(4  -  5   sin2i     -  2   cos   i)  [38] 

where  w,    the   longitude   of  perigee,   is   the   "angle"  to  ■  ft  +  w  measured  in   two 
different  planes. 

The   result  given  in  Eq.    [36]    is   sometimes   incorrectly  referred   to  as 
the   average   rate  of  rotation  of  the   line  of  apsides;   strictly  speaking  it   is 
simply   the   average   rate   of  change  of  argument   of  perigee.      The   average   rate 
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of  rotation  of  the  line  of  apsides  really  is  simply 

P-1Aeg  -  (C/2)[(2  -  3  sin2t)eQ  -  eh  sin  2t  sin  w] 
the  magnitude  of  which  is 

B*  »  (C/2)[(2  -  3  sin2t)2  +  (sin  2t  sin  u)2]1/2    [39] 

The  difference  between  this  result  and  the  value  given  in  Eq.  [36]  is  known 

to  astronomers  as  the  difference  between  the  draconitic  and  the  inertial 

motion  of  the  apsides  (1). 

It  is  of  some  interest  to  have  a  formula  for  the  mean  angular  velocity 

of  the  e  ,  e_,  e,  triad.   It  is  not  difficult  to  show  that  this  is  given  by 
s   8   h  b  j 

the  expression 

BTriad  =    (C/2)[(2   -   3  sin2t)eh  -   i    sin  2i]  [40] 

The   results   of  principal   interest  are  given  by  Equations    [26]    through 
[40] .      These  equations   or  their  equivalents   seem  to  have  been   first  obtained 
by  R.   E.    Roberson    (2)    and  by  D.   King-Hele    (3).      Both   of  them  also  discuss 
variations   in  the  orbital  period   from  the  normal  value 

P  -  2,a3/V1/2 
but   this   analysis   seems   to  be  of  a  higher  order  of  difficulty   than  that 
presented  here.     There  are  obvious  practical   difficulties   in  even  defining 
an  orbital  period.      It   could,    for  example,   refer  to  the  successive  passages 
through   the  perigee  point    (of  the   IEO  determined  at  perigee) ,   or  it   could 
refer  to  successive  passages   of  the  descending  node.      The   analysis   leading 
to  Eq.    [12]    has  been   adapted   from  a  recent  popular  textbook    (4). 

Grateful   acknowledgment  is  made  of  helpful  suggestions   from 
Dr.   W.    E.    Bleick  of  the  Naval  Postgraduate   School   and  from  Dr.   W.   W. 
King  of  Georgia  Institute   of  Technology. 
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Appendix  A. 
Remarks  about  Eq.  [3]. 

A  large  literature  has  been  developed  concerning  the  potential 
of  various  mass  distributions.   Much  of  it  has  been  concerned  with 
the  problem  of  determining  the  potential  at  an  exterior  point  of 
a  homogeneous  ellipsoid.   This  problem  is  surprisingly  difficult, 
but  if  the  figure  is  of  revolution,  the  results  can  be  expressed  in 
elementary  (even  though  quite  complicated)  terms.   Among  readily 
available  treatments,  that  of  Thomson  and  Tait  (5)  is  quite  readable 
even  though  it  is  quite  old.   In  an  editorial  footnote  added  in  1912, 
H.  Lamb  refers  to  a  history  of  the  subject  by  Todhunter.   The  latter 
seems  to  have  been  reprinted  recently,  but  the  writer  has  not  been 
able  to  find  a  copy  (6). 

However,  the  Earth  is  certainly  not   a  homogeneous  ellipsoid, 
and,  indeed,  the  finer  details  of  the  mass  distribution  of  the  Earth 
are  presently  being  revealed  by  sophisticated  analysis  of  observations 
made  on  actual  satellite  orbits.   Our  purpose  in  this  Appendix  is 
merely  to  indicate  the  genesis  of  Eq.  [3]  and  to  present  a  very 
simple  model  which  yields  Eq.  [3], 

If  0  is  the  mass  center  of  the  mass  distribution,  0  the  exterior 
point  at  which  potential  V  is  being  evaluated,  and  T  is  the  general 
location  of  mass  element  dm,  we  let  r  denote  the  distnce  OQ,  p  denote 
the  distance  OT,  a  denote  the  cosine  of  the  angle  QOT,  and  B  denote 
the  ratio  p/r.   It  is  well  known  then  that  the  potential  may  be 
expressed  in  a  series  of  integrals  (over  the  mass  distribution)  of 
Legendre  polynomials  in  a,  viz.: 

V  »  -(G/r)/IPn(a)  8n  dm 
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where  the  sum  runs  from  zero  to  infinity  and  the  symbol  P(  ) 

n 

denotes  the  Legendre  polynomial  of  degree  n.    If  the  body  is 

homogeneous  and  nearly  spherical,  with  principal  moments  of 

inertia  1,1,    and  J  ,  this  mav  be  developed  as 
x*     yy  z 

V  =  (-y/r){l+(2r2M)~1[J  +1  +1  -3(a2I  +b2I  +c2I   )]+...} 

x     y     z  x         y         z 

where  a,  h,    and  a   are  the  direction  cosines  of  line  00  with 

respect  to  the  principal  axes  (through  0).   If  we  take  the  z 

axis  as  northward,  and  consider  an  oblate  spherical  Earth  with 

semimajor  axis  A  and  semiminor  axis  B,  this  becomes 

2      2 

V  =   -  yr_1[l  +  ^f-  (1  -    3sV2)   +   ...    ] 

lOr 

and   the    coefficient   J   in   Eq .    [3]    may  be   identified  with    the   quantity 

2      2 
(A  -B   )/10.      The   deleted   terms,    indicated  by    ...    in   the   two  preceding 

equations    indicate    terms    of  higher  order   in  v      .      However,   since   the 
Earth   is   not   a  homogeneous   spheroid,    the   value   of  J   cannot  be   cal- 
culated  from  a  knowledge   of  A  =   equatorial   radius    and  B  =  polar 
radius,   but,    instead,   must  be    inferred    from  observed   perturbations 
of  satellite   orbits.      The    forgoing  has   given   an  explanation   of  the 
source   of  Eq .    [3]    without,   however,    offering   a  model   of   the  Earth 
from  which   Eq.    [3]    may  be   derived. 

Another  point   of  view  is   given   in  what    follows.      The   oblateness 
of   the   Earth  may  be    regarded   as    resulting   in   a  deficiency  of  mass 
near  the  poles    as   compared   to   a  homogeneous   spherical   Earth.      Thus, 
we   are    led   to   consider  a  model   composed   of  one  positive   point  mass 
M  located   at   0  and   two  equal  negative   point  masses    (-m)    located   on 

the  polar  axis    at  equal   distances   d  from  0;    see   Fig.    2.      The  potential 

-1  -1-1  2  2° 

V  at   0  Is   given  by  -V/G  =  My       -  mv.      -  rnr„    .      Now  v     =  x     +  z'\    so 

that  we  have 

(rl  2)?  =  x2+(z±&2  =  r-tlzd+d2  =  r2[l+r~2(d?±2zd)} 
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where   the   upper    (negative)    sign    refers    to  v     and    the    lower    (positive) 
sign   refers    to   r„.      Combining  these   relations,   we    find 

Vr/Qrr  =   -CVn)    +  2    -    (d/r)2  (1    -    3n2/r>2)   +   ... 

2 
plus    terms    involving   higher  powers    of    (d/r)    .      Since  we   propose    to 

keep    (d/r)    quite   small,   we  will   he    justified    in   neglecting    them.    Thus , 

writing  n  =  m/Mt    a   positive   number,   we   have 

V  =  -  QMr     (I  -  In)   -  nMGcTr     (1  -  3a"r     )   +   ... 

Now  M(l-2n)    =  M-2m   =  M_,  the  mass  of  the  Earth,  so  that  M   =  .V_/(l-2n). 

Thus,  recalling  that  u  =  GM_,  we  have 

h 

V  =  -  ur~l{l+  [nd2/(l-2n)]r~2(l  -   3S2r~2)   +   ...} 

2 
so   that  we   identify  J/3  with   the   quantity  nd  / (l-2n) .      Thus 

5.41xl0'4  =  JrI2/3  =   [n(l-2n)](r7/rv,)2 

and 

(dfrJ?*   5.AlxlO~4(n_1-2) 
h 

There  is,  of  course,  no  unique  decomposition.   If  we  take  n   =  1/3, 

then  m  =  M„   and  M  =  W  ,  with  ^  =  0.023*%,  =  92.3  statute  miles,  and 

E  E  E 

so  on.   The  following  tabulation  gives  a  few  such  choices. 

5 


M/ME 

m/ME 
50 

E 

J    (miles) 

101 

1.08-10"5 

13.0 

21 

10 

5.A1-10"5 

29.2 

5 

2 

2.71-10"4 

63.6 

3 

1 

5.41-10-4 

92.3 

2 

.5 

1.08-10"3 

130. A 

Since  these  values  of  (dfrS)      are  so  small,  then,  a  forteriori , 
values  of  (d/r)   will  be  even  smaller,  where  r   is  the  distance 

from  the  center  of  the  earth  to  a  satellite:  thus,  truncation 

-1/2 
of  the  binomial  series  expansions  of  (r.    _)     is  well  justi- 
fied and  any  of  these  models  does  a  good  job  of  representing 


16 


the  Earth's  actual  potential  field. 
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Appendix  B. 
List  of  Notations. 

Except  for  the  division  into  general  catagories  (general,  points, 
vectors,  scalars) ,  the  listing  is  in  the  sequence  that  the  svmhols  are 
introduced  in  the  text  rather  than  in  alphabetical  order. 

A.  General 

A  denotes  net  change  (per  orbit  passage)  of  quantity  following. 

d 
■■-   —  denotes  differentiation  with  respect  to  time,  t. 

* 

emphasizes  mean   value  over  one  orbit  passage. 

B.  "Points 

0  =  mass  center  of  the  Earth. 

0  =  position  of  satellite. 

P  =  perigee  point  of  IEO. 

D  =  descending  node  of  IEO. 

A  =  ascending  node  of  IEO. 

B  =  point  on  orbit  on  vector  1  extended. 

T  =  location  of  mass  element  dm. 

C.  Vectors 

n  =  unit  North  vector. 

a  ,D  =  fixed  unit  vectors  in  Earth's  equatorial  plane. 

r  =  vector  00  =  position  vector  of  satellite. 

V  ■  r  =  velocity  of  satellite 

'  =  attractive  force  per  unit  mass  of  satellite. 

e  =  r/r  =  unit  vector  in  direction  of  satellite  position. 
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h  =  r*V  =*  angular  momentum  about  0  per  unit  mass  of  satellite. 

e,  =  h/h   =  unit  vector  in  direction  of  h. 
h 

S  =  vector  satisfying  Eq.  [13]. 

e  =  S/s  =  unit  vector  from  0  toward  P. 
s 

j  =  e ,  *n  esc  i   =  unit  vector  from  0  toward  I). 
n 

1  =  j*e  =  unit  vector  of  triad  i,  j,  e,  . 
s  h 

ert  =  e,  xe  =  unit  vector  of  triad  e  ,  ent   e,  . 
9    h   s  s  '   6   h 

I  =  jxn  =  unit  vector  of  triad  I,  j,  n. 

B„  =  mean     angular  velocitv  of  line  of  nodes. 

N 

_* 

B„,  ,  ,  =  mean   angular  velocitv  of  e  ,  ert,  e,  triad. 
Triad  s    9   h 


D.   Scalars;  dimensionality  is  indicated  following  semicolon, 
x   =  r*a;  I. 
y   =  r»b;  L. 
z   =  r«n;  L. 

2   -2 

V  =  geopotential  at  point  Q;  £  T     . 

u  =  GM„   =  9.66x10  (statute  miles)  (seconds)   ;  L  T     . 
G  =  coefficient  of  universal  gravitation:  L^F     T     . 

J  =  coefficient  of  second  zonal  harmonic  of  potential 

-3  2   2 
=  1.623x10  V;;  L  . 


r   =  /r-r   =  distance  OQ :  L. 

-4,.    _  2  -2,   _-4 

A  =  r      (1  -  5s  r  )  ;  L   . 

-5   -A 
B  =  2zr      ;  L   . 


Also,    in  Appendix  A3    letters 
A  and  B  d,enote  semi-axes  of 
the  Earth. 


t   =  time;  7. 


2  -1 


h   =  /h»h  =  magnitude  of  angular  momentum  per  unit  mass;  r,  y   . 

y  =  true  anomaly,  angle  POQ  =  cos   (e  *e  ) ;  dimensionless . 

2  -1 
p  =  h  \i        =  semi-latus  rectum  of  IEO;  £. 

s   =  /s  «S  =  eccentricity  of  IEO;  dimensionless. 


19 
i   =  inclination  of  plane  of  TEO  with  resnert  to  equatorial  plane 

=  cos   (n*e.  ):  dimensionless. 
a)  =  argument  of  nerigee  =  angle  AOP ;  dimensionless. 
u   =  angle  H00  =  cos   (1 *e  );  dimensionless. 
w   =  angle  ROP  =  cos   (i *e  );  dimensionless. 

a  -   semimajor  axis  of  IEO;  L. 

3/2  -l/0 
P  =  2-na        u     =  period  of  orbital  passage;  T. 

1  In  o  i  In  n    -n 

=  coefficient  in  several  important  results: 

r     =  radius  of  the  Earth;  T,. 

-2  -2 

a   =  y v       =  acceleration  of  gravity  at  surface  of  Earth;  LT 

.*  -1 

u  =  mean   rate  of  change  of  argument  of  perigee;  T 

•  *  _1 

ft  =  mean    rate  of  change  of  longitude  of  ascending  node;  T 

.*  _1 

co  =  mean   rate  of  change  of  longitude  of  perigee;  T 

*  -1 

B.    -  mean   rate  of  rotation  of  line  of  apsides;  T 
A 

p  -  distance  OT  to  mass  element  dm:   L. 
a  ■  cosine  of  angle  HOT;  dimensionless. 
8  =  p/r;  dimensionless. 

P  (  )  =  Legendre  polynomial  of  degree  n;  dimensionless. 

2 
I   .1  ,1     =  principal  moments  of  inertia  at  0;  FLT   . 
X     y3    z        '      v 

a>h3c  =   direction  cosines  of  00  with  respect  to  principal  axes; 
dimensionless . 

A,B  =  semiaxes  of  oblate  Earth;  L. 

2-1 

M   =  postttve   mass  at  0;  FT  L 

2   -1 
-m   =  negat%ve   mass  on  polar  axis;  FT  L 

d   =  distance  from  M   to  -mi   L. 

n   =  m/M\    dimensionless. 

2  -1 
AC  =  M-2m   =  mass  of  Earth;  FT  L      . 
E 
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1      '        w  l 


Figure  la 


|  Figure  lb 


Figure  1   Details  of  orbit  of  Earth  satellite.   Figure  la 
is  looking  south  along  axis  of  Earth:  Figure  lb 
is  looking  along  line  of  nodes  from  ascending  node 
A  toward  descending  node  D;  Figure  lc  is  true  pro- 
jection of  orbit.   P  is  the  perigee  point,  0  is  the 
center  of  the  Earth,  and  Q  is  the  location  of  the 
satellite.  Angle  AOR  is  a  right  angle.  Q   is  the 
longitude  of  the  ascending  node,  w  is  the  argument 
of  perigee,  v   is  the  true  anomaly,  and  angles u  and  w 
are  as  indicated. 
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Figure  2   Simplified  model  of  the  F.arth 
consisting  of  positive   mass  m 
at  center  and  two  negative   masses 
-m   on  polar  axis. 

* x > 
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